*------------------------------------------------------------------------------
*Title: Political Learning, Fiscal Constraints, and the Lifetimes of Bureaus
*By:  Daniel C. Carpenter and David E. Lewis
*Date:  December 30, 2003
*
*My hope is to read the file in to stata, take care of missing data
*set the data for survival analysis, and estimate all of the models for the paper.
*------------------------------------------------------------------------------

*First, we need to do some housekeeping.
set memory 64000
set l 72
set level 95
set matsize 150

*Now, we look at the spell file.  This is based on the spell file created from the
*Lewis (2000/2003) administrative agency data set.  It selects out a few cases that are
*temporary or were found to be courts.  It adds more data to the basic file, including,
*importantly, data on surpluses adjusted for inflation.

use c:\stata\durdatacarp

*Now we describe the data.
describe

*Now we take care of some missing cases.
keep if cumdur~=.
keep if year~=.
keep if estate~=.

*We drop the temporary agencies from the analysis.

keep if temp==0

*We drop agencies that in the update process now appear as if they should have been excluded.  See update notes on Lewis website (www.princeton.edu/~delewis).

drop if agencyid==10125
drop if agencyid==10258
drop if agencyid==10213
drop if agencyid==10270
drop if agencyid==10384
drop if agencyid==10261

*We drop the courts that are left, the US Court of Appeals for the Federal Circuit and
*the US Court of Veterans Appeals

drop if agencyid==10439
drop if agencyid==10458

*We also need to change some of the agencies coded as statutorily created to executively created.  See update notes on Lewis website (www.princeton.edu/~delewis).

replace leg=0 if agencyid==10217|agencyid==10229|agencyid==10157|agencyid==10284|agencyid==10366|agencyid==10365|agencyid==10375|agencyid==10471|agencyid==10484|agencyid==10481|agencyid==10485

*Finally, we need to replace an error in the data related to the start year of the ANA.

replace startyr=. if agencyid==10554

drop y1946-y1997

tab startyr, generate (sy)

*cumdur is cumulative duration or age.  year is what is appears to be.
*estate is end state.  It is 0,1 depending on whether an agency is active
*or terminated.

*summarize
 
*First, we need to take care of one problem which is that if an 
*agency was terminated on the first of the year its cumdur is 0 instead of
*1.  We also need to take care of the problem that if an agency was
*terminated on the first of the year that its cumdur=cumdur for the 
*previous case.
replace cumdur=1 if cumdur==0
replace cumdur=cumdur+1 if cumdur[_n-1]==cumdur

*Now we try to get stata to recognize the spell file.  stset
*tells stata this is a spell file, cumdur is the agency age in days
*id is the unique agency id code.  There are multiple observations on single
*agencies.  Scale tells stata to report the results in years.

stset cumyear estate, id(agencyid) 

*Now we describe the survival-time data.

stdes
more

*Now we summarize the data

stsum
more

*First we graph the product limit estimate of the hazard rate which 
*is Figure 2.  This isn't so helpful and looks decreasing but non-monotonic.  

sts graph, hazard l1(Hazard Rate) t1(Figure 2. Non-parametric Estimates of Hazard Rate of Bureau Mortality) b2("Age of Agency in Years") b1("") title("") xtitle("")
more

*We create some of the variables that will be useful in the analyses.
*We make the difference measure an absolute value.
generate abdistnc=abs(cshmed-cspres)
generate abhsecha=abs(hsechang)
generate abpresch=abs(preschan)

*Now we generate the interaction terms.
generate partintr=unfrmaj*unfrpres*unified
generate prefinter=abhsecha*abpresch if hsecha>0&prescha>0|hsecha<0&prescha<0
replace prefinter=0 if hsecha>0&prescha<0|hsecha<0&prescha>0|hsecha==0|prescha==0
generate hunprefinter=prefinter*100

*This next version of the ideology interaction allows bigger preference divergence to decrease the importance of the other variables.
generate abdistnc2=.566-abdistnc
generate ideointr2=abhsecha*abpresch*abdistnc2

*Now interacting preference change in one direction with the equivalent of unified government, abdistnc2.
generate totalprefinter=prefinter*abdistnc2

*Now we change the party measures to Republican indicators
generate reppres=1 if dempres==0
replace reppres=0 if dempres==1
generate rephouse=1 if demh==0
replace rephouse=0 if demh==1

*This next variable is the real growth in government agencies.  Counleg is the number of legislatively created agencies.  Counpma is the
*number of agencies created by executive action.  Terminat is the number of agencies terminated in a given year.
generate realnum=counleg+counpma-terminat

*First I will estimate different models to discern which fits best.  In the first case we look at nested models.  
generate ln92budg=ln(budget)

streg realsurplusimf unified leg reorg sec line unfrmaj unfrpres partintr unemp war reppres rephouse, dist(loglog)

*To test if gamma=1 I test if ln_gam=0
test [ln_gam]_cons=0

*To test the specification we regress the dv on the linear prediction and linear prediction squared.
linktest, dist(gamma)

more
stcurve, hazard s(i) c(m) l1(Hazard Rate) t1(Log-logistic Regression) xlabel (0(5)50) b2("Age of Agency in Years") ylabel (0(.01).08) ytitle("") xtitle("") title("")        

more
streg realsurplusimf unified leg reorg sec line unfrmaj unfrpres partintr unemp war reppres rephouse, dist(gamma)

test [kappa]_cons=0
test [kappa]_cons=1
test [ln_sig]_cons=0

linktest, dist(gamma)
stcurve, hazard s(i) c(m) l1(Hazard Rate) t1(Generalized Gamma Regression) xlabel (0(5)50) b2("Age of Agency in Years")ylabel (0(.01).08) ytitle("") xtitle("") title("")

more

*Now an alternate specification just to demonstrate that in some specifications we cannot reject the null that k=0.

streg realsurplusimf abdistnc2 leg reorg sec line abhsecha abpresch totalprefinter unemp war cspres cshmed, dist(gamma)
test [kappa]_cons=0
test [kappa]_cons=1
test [ln_sig]_cons=0

stcurve, hazard s(i) c(m) l1(Hazard Rate) t1(Figure 1. Estimates of Hazard Rate from Generalized Gamma Regression) xlabel (0(5)50) b2("Age of Agency in Years")ylabel (0(.01).08) ytitle("") xtitle("") title("")

more

*Now for the AIC stuff

streg realsurplusimf unified leg reorg sec line unfrmaj unfrpres partintr unemp war reppres rephouse, dist(gamma)
predict dev, deviance
predict xb, xb
scatter dev cumyear, yline(0) m(o) t1(Generalized Gamma Residuals)
more
scatter dev xb, yline(0) m(o) t1(Generalized Gamma Residuals)
more
estimates store full
quietly streg leg reorg sec line unfrmaj unfrpres partintr unemp war reppres rephouse if xb~=., dist(gamma)
lrtest full .


streg realsurplusimf unified leg reorg sec line unfrmaj unfrpres partintr unemp war reppres rephouse, dist(gompertz)
linktest, dist(gompertz)
predict dev1, deviance
predict xb1, xb
scatter dev1 cumyear, yline(0) m(o) t1(Gompertz Residuals)
more
scatter dev1 xb1, yline(0) m(o) t1(Gompertz Residuals)
more

streg realsurplusimf unified leg reorg sec line unfrmaj unfrpres partintr unemp war reppres rephouse, dist(lognormal)
linktest, dist(lognormal)
predict dev2, deviance
predict xb2, xb
scatter dev2 cumyear, yline(0) m(o) t1(Log-Normal Residuals)
more
scatter dev2 xb2, yline(0) m(o) t1(Log-Normal Residuals)
more

streg realsurplusimf unified leg reorg sec line unfrmaj unfrpres partintr unemp war reppres rephouse, dist(loglog)
predict dev3, deviance
predict xb3, xb
scatter dev3 cumyear, yline(0) m(o) t1(Log-Logistic Residuals)
more
scatter dev3 xb3, yline(0) m(o) t1(Log-Logistic Residuals)
more
estimates store full
quietly streg leg reorg sec line unfrmaj unfrpres partintr unemp war reppres rephouse if xb~=., dist(loglog)
lrtest full .

streg realsurplusimf unified leg reorg sec line unfrmaj unfrpres partintr unemp war reppres rephouse, dist(weibull) time
linktest, dist(weibull) time
predict dev4, deviance
predict xb4, xb
scatter dev4 cumyear, yline(0) m(o) t1(Weibull Residuals)
more
scatter dev4 xb4, yline(0) m(o) t1(Weibull Residuals)
more

streg realsurplusimf unified leg reorg sec line unfrmaj unfrpres partintr unemp war reppres rephouse, dist(exponential) time
linktest, dist(exponential) time
predict dev5, deviance
predict xb5, xb
scatter dev5 cumyear, yline(0) m(o) t1(Exponential Residuals)
more
scatter dev5 xb5, yline(0) m(o) t1(Exponential Residuals)
more

stcox realsurplusimf unified leg reorg sec line unfrmaj unfrpres partintr unemp war reppres rephouse, nohr scaledsch(sca*) schoenfeld(sch*) basesurv(s) mgale(mg1)
stphtest, rank detail
predict dev6, deviance
predict xb6, xb
scatter dev6 cumyear, yline(0) m(o) t1(Cox Model Residuals)
more
scatter dev6 xb6, yline(0) m(o) t1(Cox Model Residuals)
estimates store full
quietly stcox leg reorg sec line unfrmaj unfrpres partintr unemp war reppres rephouse if xb6~=., nohr
lrtest full .

stcox realsurplusimf unified leg reorg sec line unfrmaj unfrpres partintr unemp war reppres rephouse, nohr
linktest

more
generate surptwo=3 if realsurplusimf >-24
replace surptwo=2 if realsurplusimf <=-24&realsurplusimf >-120
replace surptwo=1 if realsurplusimf <-120&realsurplusimf >-189
replace surptwo=0 if realsurplusimf <-189

drop s
drop sca*
drop sch*

stcox realsurplusimf unified leg reorg sec line unfrmaj unfrpres partintr unemp war reppres rephouse, nohr strata (surptwo) scaledsch(sca*) schoenfeld(sch*) basesurv(s)
stphtest, rank detail

drop s
drop sca*
drop sch*
drop _merge

stcox realsurplusimf unified leg reorg sec line unfrmaj unfrpres partintr unemp war reppres rephouse, nohr basesurv(s)


stcurve, survival at1(unified=1) at2(unified=0) ytitle("Cox Model Survival Probabilities") xtitle ("Age of Bureau") legend(label(1 "Unified Government") label(2 "Divided Government")) title ("") subtitle("Unified Government", position(11))
stcurve, survival at1(realsurp=0) at2(realsurp=-122) at3(realsurp=-244) ytitle("Cox Model Survival Probabilities") xtitle ("Age of Bureau") legend(label(1 "Budget Surplus (+ 1 SD)") label(2 "Mean Surplus/Deficit") label(3 "Large Deficit (- 1 SD)")) title ("") subtitle("Budget Surplus or Deficit", position(11))
*stcurve, survival at1(unified=1) at2(unified=0) ytitle("Cox Model Survival Probabilities") xtitle ("Age of Bureau") title ("") subtitle("Unified Government", position(11))
*stcurve, survival at1(realsurp=0) at2(realsurp=-122) at3(realsurp=-244) ytitle("Cox Model Survival Probabilities") xtitle ("Age of Bureau") title ("") subtitle("Budget Surplus or Deficit", position(11)) text(.55 40 "Large Deficit (-1 SD)") text(.35 43 "Mean Surplus/Deficit") text(.2 42 "Balanced Budget (+1 SD)") legend(off)

more

*Checking whether adding key variables to ancillary parameters improves model fit.  We cannot reject the null of no improvement in the log-logistic model.
streg realsurplusimf unified leg reorg sec line unfrmaj unfrpres partintr unemp war reppres rephouse , dist(loglog) ancillary(unemp realsurplusimf)
estimates store full
streg realsurplusimf unified leg reorg sec line unfrmaj unfrpres partintr unemp war reppres rephouse , dist(loglog) 
lrtest full ., force

*Checking whether adding key variables to ancillary parameters improves model fit.  We cannot (barely) reject the null of no improvement.
streg realsurplusimf unified leg reorg sec line unfrmaj unfrpres partintr unemp war reppres rephouse , dist(gamma) ancillary(unemp realsurplusimf) anc2(unemp realsurplusimf)
estimates store full
streg realsurplusimf unified leg reorg sec line unfrmaj unfrpres partintr unemp war reppres rephouse , dist(gamma) 
lrtest full ., force

*Now estimating stratified Cox models given the issues with proportionality indicated above.

generate triunemp=0 if unemp<4.45
replace triunemp=1 if unemp>=4.45&unemp<=7.45
replace triunemp=2 if unemp>7.45

stcox realsurplusimf unified leg reorg sec line unfrmaj unfrpres partintr unemp war reppres rephouse, nohr strata (triunemp) scaledsch(sca*) schoenfeld(sch*)
stphtest, rank detail

drop sca*
drop sch*

more

generate deficit=0 if realsurplusimf <-189.80
replace deficit=1 if realsurplusimf>=-189.80&realsurplusimf<=-20.48
replace deficit=2 if realsurplusimf>-20.48

stcox realsurplusimf unified leg reorg sec line unfrmaj unfrpres partintr unemp war reppres rephouse, nohr strata (deficit) scaledsch(sca*) schoenfeld(sch*)
stphtest, rank detail

drop sca*
drop sch*

*Now we test to see whether gridlock interval matters.

sort year
joinby year using gridlock

stcox realsurplusimf unified leg reorg sec line unfrmaj unfrpres partintr unemp war reppres rephouse gridlock, nohr
stcox realsurplusimf abdistnc2 leg reorg sec line abhsecha abpresch totalprefinter unemp war cspres cshmed gridlock, nohr

stcox realsurplusimf unified leg reorg sec line unfrmaj unfrpres partintr unemp war reppres rephouse grid2_n, nohr
stcox realsurplusimf abdistnc2 leg reorg sec line abhsecha abpresch totalprefinter unemp war cspres cshmed grid2_n, nohr

stcox realsurplusimf leg reorg sec line unfrmaj unfrpres partintr unemp war reppres rephouse gridlock, nohr
stcox realsurplusimf leg reorg sec line abhsecha abpresch totalprefinter unemp war cspres cshmed gridlock, nohr

stcox realsurplusimf unified leg reorg sec line unfrmaj unfrpres partintr unemp war reppres rephouse grid2_n, nohr
stcox realsurplusimf abdistnc2 leg reorg sec line abhsecha abpresch totalprefinter unemp war cspres cshmed grid2_n, nohr

*Now we generate an indicator to graph the impact of budget surplus/deficit on survival probabilities.  A 1 graphs all 
*cases below -181 (1 SD below mean) and a 0 is all cases >0 (1SD above mean).

generate deficit1=1 if realsurplusimf<-123.21
replace deficit1=0 if realsurplusimf>=-123.21

generate deficit2=0 if realsurplusimf <-294.75
replace deficit2=1 if realsurplusimf >-294.75&realsurplusimf <-189.80
replace deficit2=2 if realsurplusimf >-189.80&realsurplusimf <-120.42
replace deficit2=3 if realsurplusimf >-120.42&realsurplusimf <-20.48
replace deficit2=4 if realsurplusimf >-20.48&realsurplusimf <11.38
replace deficit2=5 if realsurplusimf >11.38

generate deficit3=0 if realsurplusimf <0
replace deficit3=1 if realsurplusimf >0

*now to use LR Tests to look at nested models
streg realsurplusimf unified leg reorg sec line unfrmaj unfrpres partintr unemp war reppres rephouse, dist(gamma)
estimates store full
streg realsurplusimf unified leg reorg sec line unfrmaj unfrpres partintr unemp war reppres rephouse, dist(exponential)
lrtest full ., force
more
streg realsurplusimf unified leg reorg sec line unfrmaj unfrpres partintr unemp war reppres rephouse, dist(weibull)
lrtest full ., force
more
streg realsurplusimf unified leg reorg sec line unfrmaj unfrpres partintr unemp war reppres rephouse, dist(lognormal)
lrtest full ., force


*Now with the other interaction term
stcox realsurplusimf abdistnc2 leg reorg sec line abhsecha abpresch totalprefinter unemp war cspres cshmed, nohr

*With budget
stcox realsurplusimf unified leg reorg sec ln92budg unfrmaj unfrpres partintr unemp war reppres rephouse if line==1, nohr

*Now interacting unemployment and budget to see if only large agencies are targeted.

generate lineunemp=line*unemp
generate budunemp=ln92budg*unemp

stcox realsurplusimf unified leg reorg sec line unfrmaj unfrpres partintr unemp lineunemp war reppres rephouse , nohr
more
stcox realsurplusimf unified leg reorg sec ln92budg unfrmaj unfrpres partintr unemp budunemp war reppres rephouse if line==1, nohr
more

more

*Now with manual controls for duration dependence
*Now with manual controls for duration dependence
generate dur1=1 if cumyear<1
replace dur1=0 if cumyear>1

generate dur2=1 if cumyear<2& cumyear>1
replace dur2=0 if cumyear>2|cumyear<1

generate dur3=1 if cumyear<3&cumyear>2
replace dur3=0 if cumyear>3|cumyear<2

generate dur4=1 if cumyear<4&cumyear>3
replace dur4=0 if cumyear>4|cumyear<3

generate dur5=1 if cumyear<5&cumyear>4
replace dur5=0 if cumyear>5|cumyear<4

generate dur6=1 if cumyear<6&cumyear>5
replace dur6=0 if cumyear>6|cumyear<5

generate dur7=1 if cumyear<7&cumyear>6
replace dur7=0 if cumyear>7|cumyear<6

generate dur8=1 if cumyear<8&cumyear>7
replace dur8=0 if cumyear>8|cumyear<7

generate dur9=1 if cumyear <9&cumyear>8
replace dur9=0 if cumyear >9|cumyear<8

generate dur10=1 if cumyear <10&cumyear>9
replace dur10=0 if cumyear >10|cumyear<9

generate dur11=1 if cumyear<11&cumyear>10
replace dur11=0 if cumyear>11|cumyear<10

generate dur12=1 if cumyear<12& cumyear>11
replace dur12=0 if cumyear>12|cumyear<11

generate dur13=1 if cumyear<13& cumyear>12
replace dur13=0 if cumyear>13|cumyear<12

generate dur14=1 if cumyear<14&cumyear>13
replace dur14=0 if cumyear>14|cumyear<13

generate dur15=1 if cumyear<15&cumyear>14
replace dur15=0 if cumyear>15|cumyear<14

generate dur16=1 if cumyear<16&cumyear>15
replace dur16=0 if cumyear>16|cumyear<15

generate dur17=1 if cumyear<17&cumyear>16
replace dur17=0 if cumyear>17|cumyear<16

generate dur18=1 if cumyear<18&cumyear>17
replace dur18=0 if cumyear>18|cumyear<17

generate dur19=1 if cumyear <19&cumyear>18
replace dur19=0 if cumyear >19|cumyear<18

generate dur20=1 if cumyear <20&cumyear>19
replace dur20=0 if cumyear >20|cumyear<19

generate dur21=1 if cumyear<21&cumyear>20
generate dur22=1 if cumyear<22&cumyear>21
generate dur23=1 if cumyear<23&cumyear>22
generate dur24=1 if cumyear<24&cumyear>23
generate dur25=1 if cumyear<25&cumyear>24
generate dur26=1 if cumyear<26&cumyear>25
generate dur27=1 if cumyear<27&cumyear>26
generate dur28=1 if cumyear<28&cumyear>27
generate dur29=1 if cumyear<29&cumyear>28
generate dur30=1 if cumyear<30&cumyear>29
generate dur31=1 if cumyear<31&cumyear>30
generate dur32=1 if cumyear<32&cumyear>31
generate dur33=1 if cumyear<33&cumyear>32
generate dur34=1 if cumyear<34&cumyear>33

replace dur21=0 if cumyear >21|cumyear<20
replace dur22=0 if cumyear >22|cumyear<21
replace dur23=0 if cumyear >23|cumyear<22
replace dur24=0 if cumyear >24|cumyear<23
replace dur25=0 if cumyear >25|cumyear<24
replace dur26=0 if cumyear >26|cumyear<25
replace dur27=0 if cumyear >27|cumyear<26
replace dur28=0 if cumyear >28|cumyear<27
replace dur29=0 if cumyear >29|cumyear<28
replace dur30=0 if cumyear >30|cumyear<29
replace dur31=0 if cumyear >31|cumyear<30
replace dur32=0 if cumyear >32|cumyear<31
replace dur33=0 if cumyear >33|cumyear<32
replace dur34=0 if cumyear >34|cumyear<33
logit estate realsurplusimf unified leg reorg sec line unfrmaj unfrpres partintr unemp war reppres rephouse dur1 dur2 dur3 dur4 dur5 dur6 dur7 dur8 dur9 dur10 dur11 dur12 dur13 dur14 dur15 dur16 dur17 dur18 dur19 dur20 dur21 dur22 dur23 dur24 dur25 dur26 dur27 dur28 dur29 dur30 dur31 dur32 dur33 dur34
estimates store full
predict p1
logit estate realsurplusimf unified leg reorg sec line unfrmaj unfrpres partintr unemp war reppres rephouse dur1 dur2 dur3 dur4 dur5 dur6 dur7 dur8 dur9 dur10 dur11 dur12 dur13 dur14 dur15 dur16 dur17 dur18 dur19 dur20 if p1~=.
lrtest full .

logit estate realsurplusimf unified leg reorg sec line unfrmaj unfrpres partintr unemp war reppres rephouse dur1 dur2 dur3 dur4 dur5 dur6 dur7 dur8 dur9 dur10 dur11 dur12 dur13 dur14 dur15 dur16 dur17 dur18 dur19 dur20

test dur1=dur5
test dur2=dur5
test dur3=dur5
test dur4=dur5
test dur6=dur5
test dur7=dur5
test dur8=dur5
test dur9=dur5
test dur10=dur5
test dur11=dur5
test dur12=dur5
test dur13=dur5
test dur14=dur5
test dur15=dur5
test dur16=dur5
test dur17=dur5
test dur18=dur5
test dur19=dur5
test dur20=dur5


more

*Now with the number of new agencies included.
generate countot=counleg+counpma
streg realsurplusimf unified leg reorg sec line unfrmaj unfrpres partintr unemp war reppres rephouse countot, dist(loglog)
streg realsurplusimf unified leg reorg sec line unfrmaj unfrpres partintr unemp war reppres rephouse countot, dist(gamma)
stcox realsurplusimf unified leg reorg sec line unfrmaj unfrpres partintr unemp war reppres rephouse countot, nohr

*Now with another real surplus var
generate urbcpiuse=urbancpi/100
generate surprealcpi=surplus/urbcpiuse

streg surprealcpi unified leg reorg sec line unfrmaj unfrpres partintr unemp war reppres rephouse , dist(loglog)
stcox surprealcpi unified leg reorg sec line unfrmaj unfrpres partintr unemp war reppres rephouse, nohr
logit estate surprealcpi unified leg reorg sec line unfrmaj unfrpres partintr unemp war reppres rephouse dur1 dur2 dur3 dur4 dur5 dur6 dur7 dur8 dur9 dur10 dur11 dur12 dur13 dur14 dur15 dur16 dur17 dur18 dur19 dur20
stcox surprealcpi unified leg reorg sec ln92budg unfrmaj unfrpres partintr unemp war reppres rephouse if line==1, nohr
stcox surprealcpi abdistnc2 leg reorg sec line abhsecha abpresch totalprefinter unemp war cspres cshmed, nohr
